function [C,Z_h,parm_h] = ml_price(S,K,r,T,Y,Z,m_Y,k_Y,m_Z,k_Z,m_Q,s_Y,s_Z,s_Q,r_12,r_13,r_23,gam,tau,tol)
% This function calculates call option price under the stochastic jump intensity model.
%
% 02-22-05 by Shu Yan.

% First setting the tolerance level of the computation.
if nargin<20
    TOL=1e-6;
else
    TOL=tol;
end

% Derive the matrices B and C.
Theta=[-0.5*gam*(gam-1), 0; 0, exp(-gam*log(1+m_Q)+0.5*gam*(gam-1)*s_Q^2)*(gam*(1+m_Q)-(gam-1)*exp(gam*s_Q^2))-1];
Pi=[m_Y;m_Z];
Lambda=[k_Y, 0; 0, k_Z];
Gamma=[s_Y^2, r_23*s_Y*s_Z; r_23*s_Y*s_Z, s_Z^2];

opts=odeset('RelTol',1e-6,'AbsTol',1e-7);
[TT,BC]=ode15s(@ml_ode_obj,[0 tau],[0 0 0 0 0],opts,Pi,Theta,Lambda,Gamma);
m=length(TT);
B=BC(m,1:2)';C=[BC(m,3) BC(m,4);conj(BC(m,4)) BC(m,5)];

% Transform the parameters under the objective measure to those under the risk neutral measure.
b=(1+m_Q)^(-0.5*gam)*exp(0.25*gam*(gam+1)*s_Q^2);
m_Q_h=(1+m_Q)*exp(-gam*s_Q^2)-1;
Z_h=b*Z;
s_Z_h=b*s_Z;
Pi_h=[1 0; 0 b]*(Pi+Gamma*B);
Lambda_h=[1 1/b; b 1].*(Lambda-gam*[r_23*s_Y 0; r_23*s_Z 0]+2*Gamma*C);
m_Y_h=Pi_h(1,1);m_Z_h=Pi_h(2,1);
k_YY_h=Lambda_h(1,1);k_YZ_h=Lambda_h(1,2);k_ZY_h=Lambda_h(2,1);k_ZZ_h=Lambda_h(2,2);
s_Y_h=s_Y;s_Q_h=s_Q;
Gamma_h=[s_Y_h^2, r_23*s_Y_h*s_Z_h; r_23*s_Y_h*s_Z_h, s_Z_h^2];

parm_h=[m_Y_h,k_YY_h,k_YZ_h,m_Z_h,k_ZY_h,k_ZZ_h,m_Q_h,s_Y_h,s_Z_h,s_Q_h,r_12,r_13,r_23];
% Set the integration upper (and lower) limit for the inverse Fourier transform.
%kmax=max(50,10/sqrt(V*min(tau)));                  
kmax=100;

% The formula next gives complex values b.c. we integrate from 0 to kmax.  Take the real part to get the true price. 
Integral=quadv('ml_price_integrand',0,kmax,TOL,[],S,K,r,T,Y,Z_h,k_YY_h,k_YZ_h,k_ZY_h,k_ZZ_h,...
    m_Q_h,s_Y_h,s_Z_h,s_Q_h,r_12,r_13,Pi_h,Gamma_h);

C=S-(exp(-r.*T).*K./pi.*Integral);              % Option price

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=ml_ode_obj(t,y,Pi,Theta,Lambda,Gamma)
%This function defines the ODEs satisfied by B and C.
%
% 02-22-05 by Shu Yan.

dy=zeros(5,1);
y12=[y(1);y(2)];y345=[y(3) y(4);conj(y(4)) y(5)];              

C_matrix=Theta+y345*Lambda+conj(Lambda')*(y345)+2*(y345)*Gamma*y345;

dy=[(conj(Lambda')+2*y345*Gamma)*y12+2*y345*Pi;...                  
    C_matrix(1,1);...
    C_matrix(1,2);...
    C_matrix(2,2)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H = ml_price_integrand(k,S,K,r,T,Y,Z,k_YY,k_YZ,k_ZY,k_ZZ,m_Q,s_Y,s_Z,s_Q,r_12,r_13,Pi,Gamma)
% This function defines the integrand of the Inverse Fourier Transform in the option price formula.
%
% 02-22-05 by Shu Yan.

ki=0.5;

k=k+ki*i;                                        % The imaginary part of k is 0.5

if length(Y)==1                                  % Expand Y and Z into vectors if they are not
    Y=repmat(Y,size(K));
    Z=repmat(Z,size(K));
end

U=[Y,Z]';

H=[];                                            % Initialize the integrand

% To reduce computation time, we only select a few contracts with different tau's to compute A, B, and C.
% The reason is that A, B, and C only depend on the parameters and maturity tau.
dT=diff([0;T]);
id_p=find(dT>0);                                 % The first contract for each maturity
K_length=diff([id_p;length(K)+1]);               % The length of each option series with identical maturity
TT=T(id_p);                                      % Different values of maturity in vector tau
TT=TT';                                          % Make it a row vector

opts=[];                                         % Set the tolerance level of the ODE solver
opts=odeset('RelTol',1e-6,'AbsTol',1e-7);

k_c=-0.5*(k.^2-i.*k);
for j=1:length(k)
    % Defining the remaining coefficient matrices.
    Theta=[k_c(j),0;0,i*k(j)*m_Q+exp(-i*k(j)*log(1+m_Q)+k_c(j)*s_Q^2)-1];
    Lambda=[k_YY-i*k(j)*r_12*s_Y,k_YZ;k_ZY-i*k(j)*r_13*s_Z,k_ZZ];
    H_j=zeros(length(K),1);                      % Initialize the integrand for a particular k                   
    
    % Solve for A, B, and C with initial values [0 0 0 0 0 0].
    [tt,ABC]=ode15s(@ml_ode_rkn,[0 TT],[0 0 0 0 0 0],opts,Pi,Theta,Lambda,Gamma);
   
    % Calculate the value of the integrand at different maturities.
    % The computation is a bit different depending whether there is one or more maturities.
    % Note that C is a conjuagate symmetric complex matrix rather than a real symmetric matrix.  We deal with it carefully.
    if length(TT)>1
        for m=2:length(tt)
            A=ABC(m,1);B=ABC(m,2:3);C=[ABC(m,4) ABC(m,5);conj(ABC(m,5)) ABC(m,6)];          
            ind=id_p(m-1):id_p(m-1)+K_length(m-1)-1;
            BU=(B*U(:,ind))';
            UCU=diag(U(:,ind)'*C*U(:,ind));
            H_j(ind)=1./(k(j)^2-i*k(j)).*exp(-i*k(j).*(r(ind).*TT(m-1)+log(S(ind)./K(ind)))+A+BU+UCU); 
        end
    else
        m=length(tt);
        A=ABC(m,1);B=ABC(m,2:3);C=[ABC(m,4) ABC(m,5);conj(ABC(m,5)) ABC(m,6)]; 
        BU=(B*U)';
        UCU=diag(U'*C*U);
        H_j=1./(k(j)^2-i*k(j)).*exp(-i*k(j).*(r.*TT+log(S./K))+A+BU+UCU);
    end
    
    H=[H,real(H_j)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% You need to cut/paste the following section into a new matlab file %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=ml_ode_rkn(t,y,Pi,Theta,Lambda,Gamma)
% This function defines the ODEs satisfied by A, B, and C.
% It is a bit different from "ml_ode_obj.m" which only considers B and C.
%
% 02-22-05 by Shu Yan.

dy=zeros(6,1);
y23=[y(2);y(3)];y456=[y(4) y(5);conj(y(5)) y(6)];              

C_matrix=Theta+y456*Lambda+conj(Lambda')*(y456)+2*(y456)*Gamma*y456;

dy=[conj(Pi')*y23+0.5*conj(y23')*Gamma*y23+[1 1]*(Gamma.*y456)*[1;1];...
    (conj(Lambda')+2*y456*Gamma)*y23+2*y456*Pi;...                  
    C_matrix(1,1);...
    C_matrix(1,2);...
    C_matrix(2,2)];


